home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / izoom.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  15KB  |  712 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  *    izoom- 
  19.  *        Magnify or minify a picture with or without filtering.  The 
  20.  *    filtered method is one pass, uses 2-d convolution, and is optimized 
  21.  *    by integer arithmetic and precomputation of filter coeffs.
  22.  *
  23.  *                     Paul Haeberli - 1988
  24.  */
  25. #include "math.h"
  26. #include "stdio.h"
  27. #include "assert.h"
  28. #include "izoom.h"
  29.  
  30. typedef struct filtinteg {
  31.     float rad, min, max;
  32.     float *tab;
  33. } filtinteg;
  34.  
  35. float flerp();
  36.  
  37. #define GRIDTOFLOAT(pos,n)    (((pos)+0.5)/(n))
  38. #define FLOATTOGRID(pos,n)    ((pos)*(n))
  39. #define SHIFT             12
  40. #define ONE             (1<<SHIFT)
  41. #define EPSILON            0.0001
  42. #define FILTERRAD        (blurfactor*shape->rad)
  43. #define FILTTABSIZE        250
  44.  
  45. static makexmap();
  46. static setintrow();
  47. static xscalebuf();
  48. static addrow();
  49. static divrow();
  50. static FILTER *makefilt();
  51. static freefilt();
  52. static applyxfilt();
  53. float filterinteg();
  54. static mitchellinit();
  55.  
  56. float filt_box();
  57. float filt_triangle();
  58. float filt_quadratic();
  59. float filt_mitchell();
  60. float filt_gaussian();
  61. float filt_lanczos2();
  62. float filt_lanczos3();
  63.  
  64. static int (*xfiltfunc)(); 
  65. static float blurfactor;
  66. int izoomdebug;
  67.  
  68. static filtinteg *shapeBOX;
  69. static filtinteg *shapeTRIANGLE;
  70. static filtinteg *shapeQUADRATIC;
  71. static filtinteg *shapeMITCHELL;
  72. static filtinteg *shapeGAUSSIAN;
  73. static filtinteg *shapeLANCZOS2;
  74. static filtinteg *shapeLANCZOS3;
  75. static filtinteg *shape;
  76.  
  77. static filtinteg *integrate(filtfunc,rad)
  78. float (*filtfunc)();
  79. float rad;
  80. {
  81.     int i;
  82.     float del, x, min, max;
  83.     double tot;
  84.     filtinteg *filt;
  85.  
  86.     min = -rad;
  87.     max =  rad;
  88.     del = 2*rad;
  89.     tot = 0.0;
  90.     filt = (filtinteg *)mymalloc(sizeof(filtinteg));
  91.     filt->rad = rad;
  92.     filt->min = min;
  93.     filt->max = max;
  94.     filt->tab = (float *)mymalloc(FILTTABSIZE*sizeof(float));
  95.     for(i=0; i<FILTTABSIZE; i++) {
  96.     x = min+(del*i/(FILTTABSIZE-1.0));
  97.     tot = tot+filtfunc(x);
  98.     filt->tab[i] = tot;
  99.     }
  100.     for(i=0; i<FILTTABSIZE; i++) 
  101.     filt->tab[i] /= tot;
  102.     return filt;
  103. }
  104.  
  105. float filterinteg(bmin,bmax,blurf)
  106. float bmin, bmax, blurf;
  107. {
  108.     int i1, i2;
  109.     float f1, f2;
  110.     float *tab;
  111.     float mult;
  112.  
  113.     bmin /= blurf;
  114.     bmax /= blurf;
  115.     tab = shape->tab;
  116.     mult = (FILTTABSIZE-1.0)/(2.0*shape->rad);
  117.  
  118.     f1 = ((bmin-shape->min)*mult);
  119.     i1 = ffloor(f1);
  120.     f1 = f1-i1;
  121.     if(i1<0)
  122.     f1 = 0.0;
  123.     else if(i1>=(FILTTABSIZE-1))
  124.     f1 = 1.0;
  125.     else
  126.     f1 = flerp(tab[i1],tab[i1+1],f1);
  127.  
  128.     f2 = ((bmax-shape->min)*mult);
  129.     i2 = ffloor(f2);
  130.     f2 = f2-i2;
  131.     if(i2<0)
  132.     f2 = 0.0;
  133.     else if(i2>=(FILTTABSIZE-1))
  134.     f2 = 1.0;
  135.     else
  136.     f2 = flerp(tab[i2],tab[i2+1],f2);
  137.     return f2-f1;
  138. }
  139.  
  140. setfiltertype(filttype)
  141. int filttype;
  142. {
  143.     switch(filttype) {
  144.      case IMPULSE:
  145.          shape = 0;
  146.          break;
  147.      case BOX:
  148.          if(!shapeBOX) 
  149.         shapeBOX = integrate(filt_box,0.5);
  150.          shape = shapeBOX;
  151.          break;
  152.      case TRIANGLE:
  153.          if(!shapeTRIANGLE) 
  154.         shapeTRIANGLE = integrate(filt_triangle,1.0);
  155.          shape = shapeTRIANGLE;
  156.          break;
  157.      case QUADRATIC:
  158.          if(!shapeQUADRATIC) 
  159.         shapeQUADRATIC = integrate(filt_quadratic,1.5);
  160.          shape = shapeQUADRATIC;
  161.          break;
  162.      case MITCHELL:
  163.          if(!shapeMITCHELL) 
  164.         shapeMITCHELL = integrate(filt_mitchell,2.0);
  165.          shape = shapeMITCHELL;
  166.          break;
  167.      case GAUSSIAN:
  168.          if(!shapeGAUSSIAN) 
  169.         shapeGAUSSIAN = integrate(filt_gaussian,1.5);
  170.          shape = shapeGAUSSIAN;
  171.          break;
  172.      case LANCZOS2:
  173.          if(!shapeLANCZOS2) 
  174.         shapeLANCZOS2 = integrate(filt_lanczos2,2.0);
  175.          shape = shapeLANCZOS2;
  176.          break;
  177.      case LANCZOS3:
  178.          if(!shapeLANCZOS3) 
  179.         shapeLANCZOS3 = integrate(filt_lanczos3,3.0);
  180.          shape = shapeLANCZOS3;
  181.          break;
  182.     }
  183. }
  184.  
  185. copyimage(getfunc,putfunc,nx,ny)
  186. int (*getfunc)(), (*putfunc)();
  187. int nx ,ny;
  188. {
  189.     int y;
  190.     short *abuf;
  191.  
  192.     abuf = (short *)mymalloc(nx*sizeof(short));
  193.     for(y=0; y<ny; y++) {
  194.     getfunc(abuf,y);
  195.     putfunc(abuf,y);
  196.     }
  197.     free(abuf);
  198. }
  199.  
  200. /*
  201.  *    general zoom follows
  202.  *
  203.  */
  204. zoom *newzoom(getfunc,anx,any,bnx,bny,filttype,blur)
  205. int (*getfunc)();
  206. int anx,any,bnx,bny;
  207. int filttype;
  208. float blur;
  209. {
  210.     zoom *z;
  211.     int i;
  212.  
  213.     setfiltertype(filttype);
  214.     z = (zoom *)mymalloc(sizeof(zoom));
  215.     z->getfunc = getfunc;
  216.     z->abuf = (short *)mymalloc(anx*sizeof(short));
  217.     z->bbuf = (short *)mymalloc(bnx*sizeof(short));
  218.     z->anx = anx;
  219.     z->any = any;
  220.     z->bnx = bnx;
  221.     z->bny = bny;
  222.     z->curay = -1;
  223.     z->y = 0;
  224.     z->type = filttype;
  225.     if(filttype == IMPULSE) {
  226.     if(z->anx != z->bnx) {
  227.         z->xmap = (short **)mymalloc(z->bnx*sizeof(short *));
  228.         makexmap(z->abuf,z->xmap,z->anx,z->bnx);
  229.     }
  230.     } else {
  231.     blurfactor = blur;
  232.     if(filttype == MITCHELL || filttype == LANCZOS2 || filttype == LANCZOS3) 
  233.         z->clamp = 1;
  234.     else
  235.         z->clamp = 0;
  236.     z->tbuf = (short *)mymalloc(bnx*sizeof(short));
  237.     z->xfilt = makefilt(z->abuf,anx,bnx,&z->nrows);
  238.     z->yfilt = makefilt(0,any,bny,&z->nrows);
  239.     z->filtrows = (short **)mymalloc(z->nrows * sizeof(short *));
  240.     for(i=0; i<z->nrows; i++) 
  241.         z->filtrows[i] = (short *)mymalloc(z->bnx*sizeof(short));
  242.     z->accrow = (int *)mymalloc(z->bnx*sizeof(int));
  243.     z->ay = 0;
  244.     }
  245.     return z;
  246. }
  247.  
  248. getzoomrow(z,buf,y)
  249. zoom *z;
  250. short *buf;
  251. int y;
  252. {
  253.     float fy;
  254.     int ay;
  255.     FILTER *f;
  256.     int i, max;
  257.     short *row;
  258.  
  259.     if(y==0) {
  260.     z->curay = -1;
  261.     z->y = 0;
  262.     z->ay = 0;
  263.     }
  264.     if(z->type == IMPULSE) {
  265.     fy = GRIDTOFLOAT(z->y,z->bny);
  266.     ay = FLOATTOGRID(fy,z->any);
  267.     if(z->anx == z->bnx) {
  268.         if(z->curay != ay) {
  269.         z->getfunc(z->abuf,ay);
  270.         z->curay = ay;
  271.         if(xfiltfunc) 
  272.             xfiltfunc(z->abuf,z->bnx);
  273.         }
  274.         bcopy(z->abuf,buf,z->bnx*sizeof(short)); 
  275.     } else {
  276.         if(z->curay != ay) {
  277.         z->getfunc(z->abuf,ay);
  278.         xscalebuf(z->xmap,z->bbuf,z->bnx);
  279.         z->curay = ay;
  280.         if(xfiltfunc)
  281.             xfiltfunc(z->bbuf,z->bnx);
  282.         }
  283.         bcopy(z->bbuf,buf,z->bnx*sizeof(short)); 
  284.     }
  285.     } else if(z->any == 1 && z->bny == 1) {
  286.         z->getfunc(z->abuf,z->ay++);
  287.         applyxfilt(z->filtrows[0],z->xfilt,z->bnx);
  288.         if(xfiltfunc)
  289.         xfiltfunc(z->filtrows[0],z->bnx);
  290.         if(z->clamp) {
  291.         clamprow(z->filtrows[0],z->tbuf,z->bnx);
  292.         bcopy(z->tbuf,buf,z->bnx*sizeof(short)); 
  293.         } else {
  294.         bcopy(z->filtrows[0],buf,z->bnx*sizeof(short)); 
  295.         }
  296.     } else {
  297.     f = z->yfilt+z->y;
  298.     max = ((int)f->dat)/sizeof(short)+(f->n-1);
  299.     while(z->ay<=max) {
  300.         z->getfunc(z->abuf,z->ay++);
  301.         row = z->filtrows[0];
  302.         for(i=0; i<(z->nrows-1); i++) 
  303.         z->filtrows[i] = z->filtrows[i+1];
  304.         z->filtrows[z->nrows-1] = row;
  305.         applyxfilt(z->filtrows[z->nrows-1],z->xfilt,z->bnx);
  306.         if(xfiltfunc)
  307.         xfiltfunc(z->filtrows[z->nrows-1],z->bnx);
  308.     }
  309.     if(f->n == 1) {
  310.         if(z->clamp) {
  311.         clamprow(z->filtrows[z->nrows-1],z->tbuf,z->bnx);
  312.         bcopy(z->tbuf,buf,z->bnx*sizeof(short)); 
  313.         } else {
  314.         bcopy(z->filtrows[z->nrows-1],buf,z->bnx*sizeof(short)); 
  315.         }
  316.     } else {
  317.         setintrow(z->accrow,f->halftotw,z->bnx);
  318.         for(i=0; i<f->n; i++) 
  319.         addrow(z->accrow, z->filtrows[i+(z->nrows-1)-(f->n-1)],
  320.                               f->w[i],z->bnx);
  321.         divrow(z->accrow,z->bbuf,f->totw,z->bnx);
  322.         if(z->clamp) {
  323.         clamprow(z->bbuf,z->tbuf,z->bnx);
  324.         bcopy(z->tbuf,buf,z->bnx*sizeof(short)); 
  325.         } else {
  326.         bcopy(z->bbuf,buf,z->bnx*sizeof(short)); 
  327.         }
  328.     }
  329.     }
  330.     z->y++;
  331. }
  332.  
  333. static setintrow(buf,val,n)
  334. int *buf;
  335. int val, n;
  336. {
  337.     while(n>=8) {
  338.     buf[0] = val;
  339.     buf[1] = val;
  340.     buf[2] = val;
  341.     buf[3] = val;
  342.     buf[4] = val;
  343.     buf[5] = val;
  344.     buf[6] = val;
  345.     buf[7] = val;
  346.     buf += 8;
  347.     n -= 8;
  348.     }
  349.     while(n--) 
  350.     *buf++ = val;
  351. }
  352.  
  353. freezoom(z)
  354. zoom *z;
  355. {
  356.     int i;
  357.  
  358.     if(z->type == IMPULSE) {
  359.     if(z->anx != z->bnx)
  360.         free(z->xmap);
  361.     } else {
  362.     freefilt(z->xfilt,z->bnx);
  363.     freefilt(z->yfilt,z->bny);
  364.     free(z->tbuf);
  365.     for(i=0; i<z->nrows; i++)
  366.         free(z->filtrows[i]);
  367.     free(z->filtrows);
  368.     free(z->accrow);
  369.     }
  370.     free(z->abuf);
  371.     free(z->bbuf);
  372.     free(z);
  373.  
  374. }
  375.  
  376. filterzoom(getfunc,putfunc,anx,any,bnx,bny,filttype,blur)
  377. int (*getfunc)(), (*putfunc)();
  378. int anx, any;
  379. int bnx, bny;
  380. int filttype;
  381. float blur;
  382. {
  383.     zoom *z;
  384.     int y;
  385.     short *buf;
  386.  
  387.     buf = (short *)mymalloc(bnx*sizeof(short));
  388.     z = newzoom(getfunc,anx,any,bnx,bny,filttype,blur);
  389.     for(y=0; y<bny; y++) {
  390.     getzoomrow(z,buf,y);
  391.     putfunc(buf,y);
  392.     }
  393.     freezoom(z);
  394.     free(buf);
  395. }
  396.  
  397. /*
  398.  *    impulse zoom utilities
  399.  *
  400.  */
  401. static makexmap(abuf,xmap,anx,bnx)
  402. short *abuf;
  403. short *xmap[];
  404. int anx, bnx;
  405. {
  406.     int x, ax;
  407.     float fx;
  408.  
  409.     for(x=0; x<bnx; x++) {
  410.        fx = GRIDTOFLOAT(x,bnx);
  411.        ax = FLOATTOGRID(fx,anx);
  412.        xmap[x] = abuf+ax;
  413.     }
  414. }
  415.  
  416. static xscalebuf(xmap,bbuf,bnx)
  417. short *xmap[];
  418. short *bbuf;
  419. int bnx;
  420. {
  421.     while(bnx>=8) {
  422.     bbuf[0] = *(xmap[0]);
  423.     bbuf[1] = *(xmap[1]);
  424.     bbuf[2] = *(xmap[2]);
  425.     bbuf[3] = *(xmap[3]);
  426.     bbuf[4] = *(xmap[4]);
  427.     bbuf[5] = *(xmap[5]);
  428.     bbuf[6] = *(xmap[6]);
  429.     bbuf[7] = *(xmap[7]);
  430.     bbuf += 8;
  431.     xmap += 8;
  432.     bnx -= 8;
  433.     }
  434.     while(bnx--) 
  435.     *bbuf++ = *(*xmap++);
  436. }
  437.  
  438. zoomxfilt(filtfunc)
  439. int (*filtfunc)(); 
  440. {
  441.     xfiltfunc = filtfunc;
  442. }
  443.  
  444. /*
  445.  *    filter zoom utilities
  446.  *
  447.  */
  448. static addrow(iptr,sptr,w,n)
  449. int *iptr;
  450. short *sptr;
  451. int w, n;
  452. {
  453.     while(n>=8) {
  454.     iptr[0] += (w*sptr[0]);
  455.     iptr[1] += (w*sptr[1]);
  456.     iptr[2] += (w*sptr[2]);
  457.     iptr[3] += (w*sptr[3]);
  458.     iptr[4] += (w*sptr[4]);
  459.     iptr[5] += (w*sptr[5]);
  460.     iptr[6] += (w*sptr[6]);
  461.     iptr[7] += (w*sptr[7]);
  462.     iptr += 8;
  463.     sptr += 8;
  464.     n -= 8;
  465.     }
  466.     while(n--) 
  467.     *iptr++ += (w * *sptr++);
  468. }
  469.  
  470. static divrow(iptr,sptr,tot,n)
  471. int *iptr;
  472. short *sptr;
  473. int tot, n;
  474. {
  475.     while(n>=8) {
  476.     sptr[0] = iptr[0]/tot;
  477.     sptr[1] = iptr[1]/tot;
  478.     sptr[2] = iptr[2]/tot;
  479.     sptr[3] = iptr[3]/tot;
  480.     sptr[4] = iptr[4]/tot;
  481.     sptr[5] = iptr[5]/tot;
  482.     sptr[6] = iptr[6]/tot;
  483.     sptr[7] = iptr[7]/tot;
  484.     sptr += 8;
  485.     iptr += 8;
  486.     n -= 8;
  487.     }
  488.     while(n--)
  489.     *sptr++ = (*iptr++)/tot;
  490. }
  491.  
  492. static FILTER *makefilt(abuf,anx,bnx,maxn)
  493. short *abuf;
  494. int anx, bnx;
  495. int *maxn;
  496. {
  497.     FILTER *f, *filter;
  498.     int x, min, max, n, pos;
  499.     float bmin, bmax, bcent, brad;
  500.     float fmin, fmax, acent, arad;
  501.     int amin, amax;
  502.     float coverscale;
  503.  
  504.     if(izoomdebug)
  505.      fprintf(stderr,"makefilt\n");
  506.     f = filter = (FILTER *)mymalloc(bnx*sizeof(FILTER));
  507.     *maxn = 0;
  508.     if(bnx<anx) {
  509.     coverscale = ((float)anx/bnx*ONE)/2.0;
  510.     brad = FILTERRAD/bnx;
  511.     for(x=0; x<bnx; x++) {
  512.         bcent = ((float)x+0.5)/bnx;
  513.         amin = ffloor((bcent-brad)*anx+EPSILON);
  514.         amax = ffloor((bcent+brad)*anx-EPSILON);
  515.         if(amin<0)
  516.            amin = 0;
  517.         if(amax>=anx)
  518.            amax = anx-1;
  519.         f->n = 1+amax-amin;
  520.         f->dat = abuf+amin;
  521.         f->w = (short *)mymalloc(f->n*sizeof(short));
  522.         f->totw = 0;
  523.         if(izoomdebug)
  524.          fprintf(stderr,"| ");
  525.         for(n=0; n<f->n; n++) {
  526.         bmin = bnx*((((float)amin+n)/anx)-bcent);
  527.         bmax = bnx*((((float)amin+n+1)/anx)-bcent);
  528.         f->w[n] = ffloor((coverscale*filterinteg(bmin,bmax,blurfactor))+0.5);
  529.         if(izoomdebug)
  530.              fprintf(stderr,"%d ",f->w[n]);
  531.         f->totw += f->w[n];
  532.         }
  533.         f->halftotw = f->totw/2;
  534.         if(f->n>*maxn)
  535.         *maxn = f->n;
  536.         f++;
  537.        }
  538.     } else {
  539.     coverscale = ((float)bnx/anx*ONE)/2.0;
  540.     arad = FILTERRAD/anx;
  541.     for(x=0; x<bnx; x++) {
  542.         bmin = ((float)x)/bnx;
  543.         bmax = ((float)x+1.0)/bnx;
  544.         amin = ffloor((bmin-arad)*anx+(0.5+EPSILON));
  545.         amax = ffloor((bmax+arad)*anx-(0.5+EPSILON));
  546.         if(amin<0)
  547.            amin = 0;
  548.         if(amax>=anx)
  549.            amax = anx-1;
  550.         f->n = 1+amax-amin;
  551.         f->dat = abuf+amin;
  552.         f->w = (short *)mymalloc(f->n*sizeof(short));
  553.         f->totw = 0;
  554.         if(izoomdebug)
  555.          fprintf(stderr,"| ");
  556.         for(n=0; n<f->n; n++) {
  557.         acent = (amin+n+0.5)/anx;
  558.         fmin = anx*(bmin-acent);
  559.         fmax = anx*(bmax-acent);
  560.         f->w[n] = ffloor((coverscale*filterinteg(fmin,fmax,blurfactor))+0.5);
  561.         if(izoomdebug)
  562.              fprintf(stderr,"%d ",f->w[n]);
  563.         f->totw += f->w[n];
  564.         }
  565.         f->halftotw = f->totw/2;
  566.         if(f->n>*maxn)
  567.         *maxn = f->n;
  568.         f++;
  569.     }
  570.     }
  571.     if(izoomdebug)
  572.      fprintf(stderr,"|\n");
  573.     return filter;
  574. }
  575.  
  576. static freefilt(filt,n)
  577. FILTER *filt;
  578. int n;
  579. {
  580.     FILTER *f;
  581.  
  582.     f = filt;
  583.     while(n--) {
  584.     free(f->w);
  585.     f++;
  586.     }
  587.     free(filt);
  588. }
  589.  
  590. static applyxfilt(bbuf,xfilt,bnx)
  591. short *bbuf;
  592. FILTER *xfilt;
  593. int bnx;
  594. {
  595.     short *w;
  596.     short *dptr;
  597.     int n, val;
  598.  
  599.     while(bnx--) {
  600.     if((n=xfilt->n) == 1) {
  601.         *bbuf++ = *xfilt->dat;
  602.     } else {
  603.         w = xfilt->w;
  604.         dptr = xfilt->dat;
  605.         val = xfilt->halftotw;
  606.         n = xfilt->n;
  607.         while(n--) 
  608.              val += *w++ * *dptr++;
  609.         *bbuf++ = val / xfilt->totw;
  610.     }
  611.     xfilt++;
  612.     }
  613. }
  614.  
  615. /*
  616.  *    filter shape functions follow
  617.  *
  618.  */
  619. float filt_box(x)
  620. float x;
  621. {
  622.     if (x<-0.5) return 0.0;
  623.     if (x< 0.5) return 1.0;
  624.     return 0.0;
  625. }
  626.  
  627. float filt_triangle(x)
  628. float x;
  629. {
  630.     if (x<-1.0) return 0.0;
  631.     if (x< 0.0) return 1.0+x;
  632.     if (x< 1.0) return 1.0-x;
  633.     return 0.0;
  634. }
  635.  
  636. float filt_quadratic(x)
  637. float x;
  638. {
  639.     if (x<-1.5) return 0.0;
  640.     if (x<-0.5) return 0.5*(x+1.5)*(x+1.5);
  641.     if (x< 0.5) return 0.75-(x*x);
  642.     if (x< 1.5) return 0.5*(x-1.5)*(x-1.5);
  643.     return 0.0;
  644. }
  645.  
  646. static float p0, p2, p3, q0, q1, q2, q3;
  647.  
  648. /*
  649.  * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  650.  * SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.
  651.  *
  652.  */
  653. float filt_mitchell(x)    /* Mitchell & Netravali's two-param cubic */
  654. float x;
  655. {
  656.     static int mitfirsted;
  657.  
  658.     if(!mitfirsted) {
  659.     mitchellinit(1.0/3.0,1.0/3.0);
  660.     mitfirsted = 1;
  661.     }
  662.     if (x<-2.0) return 0.0;
  663.     if (x<-1.0) return (q0-x*(q1-x*(q2-x*q3)));
  664.     if (x< 0.0) return (p0+x*x*(p2-x*p3));
  665.     if (x< 1.0) return (p0+x*x*(p2+x*p3));
  666.     if (x< 2.0) return (q0+x*(q1+x*(q2+x*q3)));
  667.     return 0.0;
  668. }
  669.  
  670. static mitchellinit(b,c)
  671. float b, c;
  672. {
  673.     p0 = (  6.0 -  2.0*b         ) / 6.0;
  674.     p2 = (-18.0 + 12.0*b +  6.0*c) / 6.0;
  675.     p3 = ( 12.0 -  9.0*b -  6.0*c) / 6.0;
  676.     q0 = (       8.0*b + 24.0*c) / 6.0;
  677.     q1 = (      - 12.0*b - 48.0*c) / 6.0;
  678.     q2 = (         6.0*b + 30.0*c) / 6.0;
  679.     q3 = (       -     b -  6.0*c) / 6.0;
  680. }
  681.  
  682. #define NARROWNESS    1.5
  683.  
  684. float filt_gaussian(x)
  685. float x;
  686. {
  687.     x = x*NARROWNESS;
  688.     return (1.0/exp(x*x) - 1.0/exp(1.5*NARROWNESS*1.5*NARROWNESS));
  689. }
  690.  
  691. float filt_lanczos2(x)
  692. float x;
  693. {
  694.     if (x == 0)
  695.     return 1;
  696.     if ((x = fabs(x)) >= 2)
  697.     return 0;
  698.     x *= M_PI;
  699.     return sin(x) * sin(x / 2) / (x * x / 2);
  700. }
  701.  
  702. float filt_lanczos3(x)
  703. float x;
  704. {
  705.     if (x == 0)
  706.     return 1;
  707.     if ((x = fabs(x)) >= 3)
  708.     return 0;
  709.     x *= M_PI;
  710.     return sin(x) * sin(x / 3) / (x * x / 3);
  711. }
  712.